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Abstract 

The dynamics of magnetic reversal process plays an important role in the design of 
the magnetic recording devices in the long time scale limit. In addition to long time 
scale, microscopic effects such as the entropic effect become important in magnetic 
nano-scale systems. Many advanced simulation methods have been developed, but 
few have the ability to simulate the long time scale limit and to accurately model 
the microscopic effects of nano-scale systems at the same time. We develop a new 
Monte Carlo method for calculating the dynamics of magnetic reversal at arbitrary 
long time. For example, actual calculations were performed up to lO^'^ Monte Carlo 
steps. This method is based on microscopic interactions of many constituents and 
the master equation for magnetic probability distribution function is solved sym- 
bolically. 



1 Introduction 



Studies on complex systems with meta-stable states are valuable in many 
disciplines. For instance, aging phenomena affect the properties of materials, 
magnetic relaxation processes determine the functions of magnetic recording 
media, and the study of first order phase transitions and spinodal decompo- 
sition have wide applications. Often, slow dynamics plays a central role in 
the understanding of these complex systems. When no analytic solutions are 
available, numerical simulations are the only tools to study these systems. In 
a Monte Carlo simulation, the meta-stable states could manifest as apparently 
stable or "equilibrium" states. Only after a very long time the simulation could 
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overcome the free energy barrier and go into equilibrium. When using the stan- 
dard Metropohs algorithm, this time scale may become unattainable even with 
the fastest computer available. Hence advanced Monte Carlo techniques such 
as the Absorbing Markov Chains algorithm [1], which is a generahzation of 
the N-fold way method [2], and the Projection Method [3] have been devel- 
oped to address this issue. These methods have been used to study magnetic 
reversal [4,5,6]. In some cases, where the difficulties of mcta- stability are not 
too serious, conventional Monte Carlo methods are also used [7,8,9,10]. 

In this paper, we will present a new method of calculating the dynamics of 
systems with meta-stable states by solving the master equation for magnetic 
probability distribution function symbolically. We illustrate our method in the 
context of magnetic reversal on the two dimensional Ising square lattice. 



2 Method 



The Ising model is a basic model to study the phase transition, and has also 
been used as a standard benchmark model for testing new algorithms 
[4,11,12,13,14,15,16,17,18,19,20]. The Hamiltonian for the Ising model is 

n = -J SiSj-hY,Si = E-hM (1) 

<i,j> i 



where Si — ±1. The first summation is over nearest neighbors; the second 
summation is over all lattice sites. J and h are the exchange constant and 
external field respectively. We define two symbols E and M as the energy of 
the exchange part and the magnetization part. 



SiS 



J 



<i,j> 



M = Y.s, (2) 

i 

Ideally, to study the dynamics of magnetic reversal on the Ising model, the 
microscopic master equation, 

= Y.^{a\a')P{a\t)-u{a'\a)P{a,t) (3) 



is to be solved, a denotes the microscopic state of the system, uj{a\a') is the 
transition rate from a' to cr, and P{<J,t) is the probability distribution of the 
state a at time t. For a lattice with N sites, Eq. (3) is a set of 2^ coupled 
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Fig. 1. Log of probability distributions plotted at t = 10, 1000, 10^" and 10^^ for 
L = 10,T = O.UTc and h = 0. Initial condition is P{t = 0,M) = Sm -wo- At time 
t = 10^^, the probability distribution is almost indistinguishable from the equilib- 
rium distribution. The little spikes marked with arrows are due to band structure of 
the Ising model (Fig. 2). Equilibrium probability distribution was generated using 
the Wang-Landau method up to correction factor of log(/) = 1.25 x 10~^. Error 
bars are smaller than the thickness of the lines. 



differential equations and solving them is infeasible. Instead the dynamics 
can be approximated by writing the master equation in terms of the total 
magnetization M, 

"^^^j^'^^ = J^u{M\M')P{M', t) - u{M'\M)P{M, t) (4) 



which is a set of + 1 differential equations. The approximate form of the 
macroscopic transition rates uj{M\M') was discussed by Lee et.al. [21]. We 
made identical approximations in the transition rates as the "mean field dy- 
namics 2" (MFD2) obtained by Lee et.al. [21]. Other ways of obtaining the 
transition rates, such as the Transition Matrix Monte Carlo method [22], may 
also be used. Our main contribution in this paper is in the way of solving the 
macroscopic master equation (Eq. (4)). Our results are accurate to within the 
validity of basic approximations in the transition rates. However, the gain in 
efficiencies and accuracies of our results are many orders of magnitude larger 
than most Monte Carlo methods. 

As in Ref. [21] the equilibrium distribution is used to calculate the transition 
rates. With this approximation, not all microscopic information is available. 
The region of validity for this approximation will be discussed later. 
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An efficient Monte Carlo algorithm to calculate the density of states has been 
recently proposed by Wang and Landau [19,20]. We have used the Wang- 
Landau [19,20] algorithm to estimate the density of states as a function of 
exchange energy and magnetization g{E,M), with E and M defined by Eq. 
(2). For a given /3 = l/k^T and h, the joint probability distribution P^'^f^{E, M) 
is then obtained using Boltzmann weights as 



P;%{E, M) = g{E, M) exp(-/5[E - /iM])/Z(/3, h) 



(5) 



Z is the partition function given by, 

h)=Y. g{E, M) exp(-/3[£; - hM]) 



(6) 



EM 



Finally, the equilibrium probability distribution of M can be obtained, 

Pl%{M)^Y.n\iE,M) (7) 



E 



For a lattice with E^ = N sites, there are + 1 possible values of magnetiza- 
tion. Then the master equation (Eq. (4)) can be written in a matrix form. 



dP{t) 
dt 



= A ■ P{t) 



(8) 



where the {N -\- 1) dimensional vector P{t) is 



^ P{M ^ -N,t)^ 



P{M^O,t) 



P{M = +N,t) 



(9) 



In this paper, we present a way to calculate the explicit solution to Eq. (8). 
We restrict the transitions to | AM| = 2. This is a vital component of our work 
and we shall elaborate its implications and the conditions in which it is valid. 
Transitions which involve |AM| = 2 are all single spin fiip moves or any move 
which flip spins from +1 to —1 and A^ + 1 spins from —1 to +1 (or vice 
versa). At low temperature and low field, we do not expect to observe many 
cluster fiips or multiple spin fiip transitions. Hence our assumption is valid at 
low temperature and low field. Incidentally, the process of magnetic reversal at 
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Fig. 2. The upper figure shows a band structure that is the origin of spikes observed 
in the probabihty distribution of Fig. 1. These are local energy minimum config- 
urations, and any single spin flip increases the energy. The lower figure shows the 
domain of accessible state for L = 10 Ising model. The minimum energy state with 
M = is E = —160 or E = Eg + 4L, where Eg is the ground state energy. Arrows 
indicate band structure states shown in the upper figure, and their exact density of 
states can also be trivially calculated. 

low temperature and field is via the single-droplet process [23,24,25,26]. This 
is also the domain where the approximation to our transition rates is valid. At 
high temperatures, other magnetic reversal processes such as the multi-droplet 
process dominate and our approximation becomes not accurate in this regime. 
Short time scale behavior becomes important at high temperatures and other 
Monte Carlo methods may be used. 

With the constraint |AM| = 2, the transition rate uj is written as 



uj{Mj\Mi) 



P'"^{Mj)/[P^%Mj) + P^^Mi)] \i-j\<l 



(10) 







otherwise 
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if the Glauber rate is employed. The Metropolis rate can also be used, and we 
compare the results obtained by both rates in Table 1. 

For the transition rates given by Eq. (10), the master equation (Eq. (8)) can 
be reduced into one linear differential equation with constant coefficients and 
the analytic solutions for this equation is known [27]. A procedure to derive 
the solution is given in an appendix. The general solution is of the form, 

N 

P{t) = Y,aiVieMKt) (11) 

i=0 



where and Vi are the eigenvalues and eigenvectors of the matrix A and ctj 
are constant coefficients determined by solving Eq. (11) with initial conditions. 



N 



P{Q)=Y.aSi (12) 



i=0 



We used MATHEMATICA to calculate ctj, Vi and Aj numerically, and P{t) is 
represented symbolically as a function of time. We have also made use of a 
MATHEMATICA function for setting arbitrary precisions and compared our 
results for calculations with different precisions. With P{t) evaluated symbol- 
ically, other magnetic quantities, such as M{t), can be evaluated symbolically 
as well. 

The initial conditions were set with all the spins in the —1 position and the 
external field always favors spin configurations with Sj = +1. Our simulations 
were performed on a dual 2 GHz PowerMac G5 PC, with over 95% of the 
computing resources devoted to calculating the equilibrium density of states 
using the Wang-Landau method [19,20]. 



3 Results 

All eigenvalues are real and negative except for the largest eigenvalue which 
is zero. The probability distribution can be re-written as, 

N 

P{t)^aoVo + ^aiViex.-p{-\Xi\t) (13) 

i=l 

N 

^P"^ + J2a,VieM-\Mt) (14) 

i=l 
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Fig. 3. Magnetization versus time curve for L = 10, T = CllTc, h/J = 0.25. Equi- 
librium probability distribution was estimated using the Wang-Landau method up 
to a correction factor log(/) = 1.25 x 10""^. In agreement with previous works [25,26], 
the switching duration is approximately equal to the lifetime. Error bars are much 
smaller than thickness of the line. 




J/T 

Fig. 4. Plot of switching time r at h/ J = 0.75, which shows very good agreement 
with previous studies by Novotny [4]. Lattice sizes are L = 10 (circles) and L = 50 
(squares). Triangles show the results reported by Novotny [4] for L = 10 lattice. 



where ckq and Vq are the coefficients and eigenvector corresponding to the 
largest eigenvalue (which is zero), and a^VQ — P^^. Hence, 



lim P{t) = (15) 

t— >oo 
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Fig. 1 shows the time evohition of the probabihty distribution function for 
L = 10, T = 0.44T,, and h = 0. The magnetization was set to —N at t = and 
allowed to evolve via the master equation. At t = 10, probability distribution 
falls off exponentially from near 1 at M = -100 to exp(-225) at M = 100. At 
t — 10^°, the probability distribution at M = 100 increases and we begin to see 
two peaks, but the peak at M = 100 is 20 orders of magnitude smaller. At t — 
10^^, the probability distribution is almost symmetric and indistinguishable 
from the equilibrium probability distribution. Spikes indicated with arrows 
originate from special band structures of the Ising model (Fig. 2). These band 
structures are local energy minimums, any single spin flip from these structures 
resulted in an increase of energy. Fig. 2 also shows the domain of accessible 
states for L = 10, and arrows indicate band structures states. They have 
equal spacing in magnetization of AM = 2L and their energies are Eg + 4L, 
where Eg is the ground state energy. We have demonstrated our calculation 
of the probability distribution with zero external field, because the spikes in 
the probability distribution appear most prominently with h — 0. 

We are also interested in studying magnetic switching in an external field 
where a direct comparison with previous work is available. Fig. 3 shows reversal 
of magnetization in an external field at T = O.llTc ^-nd h/J — 0.25 for L — 
10 lattice. In agreement with previous works [25,26] the average switching 
duration is approximately equal to the lifetime. 

The temperature dependence of the switching time can be studied by plotting 
the switching times versus the inverse temperature at h/J — 0.75 (Fig. 4). 
Here we define the switching time as the time when the average magnetization 
reaches zero ((M(t)) = 0). The data for L = 10 are compared with those by 
Novotny [4]. We obtained very good agreement with previous works. In the 
insert of Fig. 4, we plot the switching time x system size as a function of 
the inverse temperature. The data for L = 10 and L = 50 collapse into a 
single curve, as expected in the single-droplet regime. In contrast with the 
other method employed in Ref. [3] where r is calculated directly, the present 
method calculates the entire probability distribution of the master equation 
that allows us to calculate other quantities such as average magnetization. 

There are two sources of approximation in our method; the first one is the 
approximation of transition rates using equilibrium distribution, the second 
one originated from uncertainty in the density of states estimated from Wang- 
Landau method. The effects of uncertainties on a L = 10 lattice in a field 
ol h/J — 0.75 are summarized in Table 1. In the Wang-Landau algorithm, 
the correction factor / is used in the process of refining the density of states. 
The final value of / is related to the accuracy of the calculation, and smaller 
"log(/)" values corresponds to more accurate equilibrium distributions. Our 
results suggest that running the Wang-Landau algorithm to log(/) value of 
1.25 X 10"^ is sufficient to within our simulation precision. Table 1 also shows 
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Table 1 

Comparing results using Glauber and Metropolis transition rates and using different 
estimates of equilibrium distributions. "log(/) values" given in the third column are 
the values of correction factors in Wang-Landau method [19] on L = 10 lattice with 
h/J = 0.75. 

the dependence of switching time r with Glauber and Metropohs rates. The 
Metropolis rates are given by, 



Lj{M\M') a; min 

There is no significant difference in the switching time between Glauber and 
Metropolis transition rates for low temperature (J/T = 2.67). But they differ 
at high temperature because the approximation to the transition rates be- 
comes less accurate at high temperatures. The uncertainties in the switching 
time, At was obtained from several independent runs. 



for |i — j| < 1 



(16) 



4 Conclusion 



In conclusion, we derived an explicit expression for the magnetic probabil- 
ity distribution as a function of time. With the probability distribution, all 
dynamical information is available. And we used the information to obtain 
consistent results with previous works. Another advantage is that dynami- 



9 



cal properties can be evaluated to arbitrary long time without requiring more 
computational resources. Lastly, we should mention that our method is general 
and not restricted to discrete systems or any specific models. Future develop- 
ment of the present method should focus on calculation of larger lattice sizes. 
Currently, lattice sizes are limited by the efficiency of Wang-Landau algorithm 
to calculate joint density of states. 
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Appendix : Solutions to the Master Equation 



A general procedure to solve the master equation will be presented in this ap- 
pendix. Firstly, we shall introduce the mathematical notation, the probability 
distribution function P{M, t) will be denoted in a form of a vector. 



P{M ^ -N + 2,t) 



\P{M = +N ,t) J 



(17) 



Note that for an Ising model with sites, with N being an even number, the 
number of possible magnetizations is + 1 with M = {—N, • • • 0, • • • N}. The 
master equation can be written in a matrix form. 



dP{t) 
dt 



A ■ P{t) 



where ^ is a (N + 1) x (A^ + 1) tridiagonal matrix with constant matrix 
elements. Suppose S~^AS = V, where 5 is a similarity transform and P is a 
diagonal matrix. Then we have 



dPjt) 
dt 



svs-^ ■ P{t) 



(19) 
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Defining a new vector Q{t) 



— S ^P{t), we obtain 



(20) 



Since "D is a diagonal matrix, Eq. (20) is easily solved with Qi{t) = Ci exp(Ait), 

— * — * 

where Aj is the ith eigenvalue of A. With P{t) — SQ{t), we have 

P(i) = exp(Aii) (21) 



We write P in terms of the eigenvectors so that the constant coefficients ctj 
can be easily found by initial conditions, 

P{t^O)^Y.^i'^i (22) 
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